from __future__ import division
from math import *
from Init_Cond import *

tstep = []  #viscosity stability timestep
tstep2 = [] #torque stability timestep

def f(M):
    del tstep[:]
    for a in range(int(N)):
        if nu[a] < 1.e-6:   #if nu is close to zero, stability equation blows up.  If nu is zero, tstep can be infinity
            tstep.append(10.*deltaT)    #just set it to be above deltaT
        else:
            tstep.append(deltaX**2./6.*r[a]/nu[a])  #stability value

def o(deltat):
    t1 = min(tstep)  #what's the lowest stable timestep?
    if len(Sat_tstep) > 0:
        t2 = min(Sat_tstep)
    else:
        t2 = deltaT
    # print(t1, t2)

    t = min(t1, t2)


    if t <= deltat: #if new timestep is smaller than previous, use newer
        return t*0.03   #guarantee it's small enough
    elif t > deltaT:    #if new timestep is larger than that set to run in Init_Cond, use deltaT
        return deltaT
    elif t > deltat and t <= deltaT:    #if new timestep is larger than previous, use newer (adapts to run code faster)
        return t*0.03


def sats(deltat, M):
    delta_t = deltat
    for a in range(len(Sat_M)):
        if Sat_M[a] > 0.0:
            j = int(round((Sat_r[a] - r_I)/deltar - .5))


            if Sat_r[a] < Sat_r_default[j]:
                j = j - 1

            for i in range(j):
                fprime = 2.5
                delta_x1 = sqrt(r[i] + deltar/2.)/(Sat_r[a] - r[i] - deltar)**3. - sqrt(r[i] - deltar/2.)/(Sat_r[a] - r[i])**3.
                delta_x2 = sqrt(r[5])/(Sat_r[a] - r[5] + deltar/2.)**3. - sqrt(r[4])/(Sat_r[a] - r[4] + deltar/2.)**3.
                A = 6.*r[i]**1.5*pi*(G*M)/(fprime*Sat_r[a]**7.*Sat_w[a]**2.)*(M/Sat_M[a])**2.
                delta_t1 = A*(Sat_r[a] - r[i] - deltar/2.)**4./(Sat_r[a] + 5.*r[i] - deltar/2.)*delta_x1

                if delta_t1 <= delta_t:
                    delta_t = 0.8*delta_t1

    return delta_t
